SOCR ≫ TCIU Website ≫ TCIU GitHub ≫

1 Background

In this vignette, we introduce the use of fmri_stimulus_detect function on our simulated brain fMRI data (generated by fmri_simulate_func) to get the p value to locate the stimulated parts in the brain. We also demonstrate the use of fmri_post_hoc on p-values stored in our package to do the FDR correction and spatial clustering. Finally, the fmri_2dvisual and fmri_3dviusal functions will be introduced to make the 2D and 3D plot of the fMRI data to better visualize the brain and its stimulated parts. The plots comparison function is also introduced to visualize the difference of stimulated parts provided by two different 3D array p value data.

require(TCIU)

2 Load and display the data

Here we demonstrate the 2D plane of our 3D nifti mask from each of the three direction. We can tell from the plot that it is a brain shell.

par(mfrow=c(2,2))
image(mask[,30,], asp=1)
image(mask[30,,], asp=1)
image(mask[,,30], asp=1)

As we can see, pval is a 3D array with its dimension as 64 times 64 times 40. Its range is between 0 to 1.

dim(pval1)
## [1] 64 64 40
summary(pval1)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.4514  1.0000  0.7882  1.0000  1.0000

3 Analysis of simulated fMRI data

Here we simulate our 4D fMRI data with its first three dimension as 64 times 64 times 40 and its fourth dimension as the 160-time length. We specify the dimension of the data, starting time points when the fMRI data receives the stimulus and the duration of each stimulated time periods to get it. As we want our data shaped like a brain, we provide a brain-shaped mask for it. Otherwise, it will generate a sphere mask automatically.

fmri_generate = fmri_simulate_func(dim_data = c(64, 64, 40), mask = mask, 
                                   ons = c(1, 21, 41, 61, 81, 101, 121, 141), 
                                   dur = c(10, 10, 10, 10, 10, 10, 10, 10))

# the output include simulated fMRI data, its mask, 
# the starting time points of the stimulated period and its duration 
# as well as all the stimulated time points
dim(fmri_generate$fmri_data)
## [1]  64  64  40 160

Here we apply stimulus_detect function on our simulated fMRI data to get the 3D array p value data. We use t-test here for simplicity. But we also have Wilcoxon test, HotellingT2 test, gLRT test for complex data, generalized likelihood test, generalized linear model and on-off intensity difference method to generate the p value to detect the stimulated parts in the fMRI data. The smaller the p value is, the higher level the stimulated part is.

p_simulate_t_test = 
fmri_stimulus_detect(fmridata= fmri_generate$fmri_data, 
                     mask = fmri_generate$mask,
                     stimulus_idx = fmri_generate$on_time,
                     method = "t-test" , 
                     ons = fmri_generate$ons, 
                     dur = fmri_generate$dur)

dim(p_simulate_t_test)
summary(p_simulate_t_test)

From the dimension of the generated p_simulate_t_test (pval2), we know it is a 3D array p value with its dimension as 64 times 64 times 40.

4 Analysis of real fMRI data

As we save a 3D array p value data generated from the real fMRI brain data, here we can apply post-hoc test on it to improve the performance. We can apply FDR correction and spatial clustering to our 3D array p value data to prevent the high proportion of false positive voxels and keep only spatially connected significant voxels. Notice that when doing spatial clustering, the threshold(spatial_cluster.thr) and cluster size(spatial_cluster.size) are needed to filter contiguous clusters of locations in a 3D array that are below some threshold and with some minimum size.

# do the FDR correction
pval_fdr = fmri_post_hoc(pval1 , fdr_corr = "fdr",
                         spatial_cluster.thr = NULL,
                         spatial_cluster.size = NULL, 
                         show_comparison = FALSE)

# do the spatial clustering
pval_posthoc = fmri_post_hoc(pval_fdr, fdr_corr = NULL,
                             spatial_cluster.thr = 0.05,
                             spatial_cluster.size = 5, 
                             show_comparison = FALSE)

5 Visualization and Comparison of fMRI data

5.1 2D and 3D visualization

We first generate the 2D plot and 3D interactive of our simulated fMRI data based on the p value generated above from sagittal, coronal and axial view.

for(axis in c("x", "y", "z")){
  axis_i = switch(axis, 
                  "x" = {35},
                  "y" = {30},
                  "z" = {22})
  print(fmri_2dvisual(pval2, list(axis, axis_i), 
                      hemody_data=NULL, mask=fmri_generate$mask, 
                      p_threshold = 0.05, legend_show = TRUE, 
                      method = "scale_p",
                      color_pal = "YlOrRd", multi_pranges=TRUE))
}

fmri_3dvisual(pval2, fmri_generate$mask, 
                            p_threshold = 0.05, method="scale_p", multi_pranges=TRUE)$plot

5.2 Comparison of performance of different methods on the same fMRI data

We can use fmri_pval_comparison_3d to visualize two p value data to see how they perform differently on detecting stimulated parts by 3D plot. Here we compare the difference of stimulated parts of two different fMRI data with the same mask, since our original fMRI is too big to use here for different statistical test.

# the two p value are the p value generated based on the simulated fMRI
# and the p value saved in the package and finished post hoc test
fmri_pval_comparison_3d(list(pval2, pval_posthoc), mask, 
                                list(0.05, 0.05), list("scale_p", "scale_p"), 
                                multi_pranges=FALSE)

We can also use fmri_pval_comparison_2d to visualize whatever number of p value (generated by different statistical tests on the same fMRI data) to see their difference by 2D plot. For simplicity here we only compare two different 3D array p value data.

fmri_pval_comparison_2d(list(pval2, pval_posthoc), 
                        list('pval_simulated', 'pval_posthoc'),
                        list(list(35, 33, 22), list(40, 26, 33)), 
                        hemody_data = NULL, 
                        mask = mask, p_threshold = 0.05, 
                        legend_show = FALSE, method = 'scale_p',
                        color_pal = "YlOrRd", multi_pranges=FALSE)